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ABSTRACT 

A simple numerical scheme is presented for the construction of three-integral phase- 
space distribution functions for oblate galaxy models with a gravitational potential 
of Stackel form, and an arbitrary axisymmetric luminous density distribution. The 
intrinsic velocity moments can be obtained simultaneously with little extra effort. 
The distribution of the inner and outer turning points of the short-axis tube orbits 
that are populated can be specified freely, and is chosen in advance. The entire 
distribution function is then derived from the density by an iterative scheme that 
starts from the explicitly known distribution function of the thin-orbit (maximum 
streaming) model, in which only the tubes with equal inner and outer turning points 
are populated. The versatility and limitations of this scheme are illustrated by the 
construction of a number of self-consistent three-integral flattened isochrone models 
of Kuzmin-Kutuzov type, and by investigation of special cases where the scheme is 
tractable analytically. This includes the behaviour of the distribution functions in 
the outer regions of the models. The scheme converges rapidly for models containing 
orbits with ratios of the outer to inner turning point as large as ten, and is particularly 
suited for the construction of tangentially anisotropic flattened models, self-consistent 
as well as non-consistent. The algorithm simplifies in the disk and spherical limit, and 
can be generalized to triaxial models. 

Key words: stellar dynamics - galaxies: kinematics and dynamics - galaxies: struc- 
ture 



1 INTRODUCTION 

The observable properties of elliptical galaxies indicate that 
their internal dynamics is governed by three integrals of mo- 
tion (Binney 1976, 1978). For oblate systems two of the 
three are known, the energy E and the angular momentum 
component L z along the symmetry axis. An exact third 
integral ^3 exists only for special classes of potentials, but 
adequate approximations have been derived for moderately 
flattened axisymmetric models (e.g., Saaf 1968; Innanen & 
Papp 1977; Gerhard & Saha 1991). 

The construction of the full class of dynamical models 
for elliptical galaxies is a major undertaking. Progress has 
been made recently on a number of fronts, in particular for 
oblate systems. Even though elliptical galaxies as a class 
have triaxial shapes, the majority may well be nearly oblate 
(Franx, Illingworth & de Zeeuw 1991), so that oblate models 
are useful. Various practical methods have been developed 
for the construction of the special model with phase-space 
distribution function / = f(E, L z ) (Hunter & Qian 1993; 
Dehnen & Gerhard 1994; Magorrian 1995; Kuijken 1995; 
Qian et al. 1995). 

An exact third integral is known explicitly for the class 
of flattened models with a potential of Stackel form (Kuzmin 
1956; de Zeeuw 1985, hereafter dZ), and some self-consistent 
three-integral dynamical models of this type have been con- 
structed, e.g., by numerical methods (Bishop 1986; 1987) 
or by series expansions (Dejonghe & de Zeeuw 1988, here- 



after DZ). The distribution function for the model with the 
maximum possible streaming motions can be found by a 
single quadrature over the density (Bishop 1987; de Zeeuw 
& Hunter 1990, hereafter ZH). In oblate Stackel models all 
orbits are short-axis tubes, but only those with vanishing 
radial action — which lie on spheroidal shells — are popu- 
lated in the maximum streaming model. They are often re- 
ferred to as thin (tube) orbits, and the corresponding model 
is called the thin-orbit model. These flattened models con- 
nect the sphere made exclusively of circular orbits with the 
similar axisymmetric disk. 

When no exact ^3 is known, dynamical models can be 
constructed by numerical methods (e.g., Richstone 1980, 
1982, 1984; Levison & Richstone 1985a, b) or by use of 
an approximate integral (Petrou 1983a, b). This approach 
has been employed recently by Dehnen & Gerhard (1993), 
who constructed a large family of approximate three-integral 
distribution functions for a flattened isochrone model, and 
investigated the relation between the internal dynamics and 
the observable kinematics. Their method is applicable to a 
wide variety of mass models with realistic density profiles. 
The one application that has been published so far is for a 
mass model that is nearly identical to the Kuzmin-Kutuzov 
model. This has a Stackel potential, and its exact third in- 
tegral has been used to construct a number of distribution 
functions (DZ, ZH). 

Little is known about the stability of flattened galaxy 
models. Some N-body simulations have been carried out 
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(Merritt 1987; Merritt & Stiavelli 1990), but the paucity of 
available distribution functions to set up the initial condi- 
tions is one of the main reasons for our lack of knowledge. 
We have shown recently (Robijn & de Zeeuw 1995) that 
the linear stability analysis pioneered by Kalnajs (1977) for 
axisymmetric disks, and subsequently used by e.g., Poly- 
achenko & Shukhman (1981), Palmer & Papaloizou (1987), 
Weinberg (1989, 1991) and Saha (1991, 1992) to study 
spherical models, also can be carried out for oblate Stackel 
models. One of the first applications is a study of the thin 
orbit models, which have been shown by N-body simulations 
to be liable to ring- and lopsided instabilities, depending on 
the flattening of the model. Based on studies of spheres and 
flat disks, we expect that an increase in the amount of radial 
support will stabilize the radially 'cold' thin-orbit models. 
In order to investigate this, we need distribution functions 
for models in which not only the thin short-axis tubes are 
populated, but also 'thick' tube orbits with a finite radial 
extent. It is those models that we construct here. 

The thin-orbit model has a distribution function of the 
form / = ftsm( Jcf>, Jv)6( J\), where J\ is the radial action, 
Jcf> = L z is the azimuthal action, and J v is the latitudinal 
action. In this paper we write the distribution function in 
the (general) form / = / gS m(J^, Jv)g(J\, Jcf>, Jv), where g 
is a preassigned function, and we show how to find / gS m, 
consistent with a given axisymmetric density p in an oblate 
Stackel potential V, by an iterative method, starting with 



the thin-orbit function /tsm as a first guess for f g: 



We 



will consider functions g that are peaked in J\, so that the 
models will be fairly close to the thin-orbit model, and few 
iterations are needed. The stability analysis of these models 
will be discussed in a subsequent paper. 

Our method of specifying part of the distribution func- 
tion, and solving for the remainder, is not new, and was 
used for flat disks by Shu (1969). Bishop (1986) applied 
it to oblate Stackel models, starting from a different initial 
guess. Gerhard (1991) and Dehnen & Gerhard (1993) have 
recently popularized this approach for spherical and oblate 
models. 

In Section 2 we define our notation, and present the con- 
struction method. A detailed description is given in Section 
3, where we also investigate what properties of the assigned 
function g are important for convergence of the iterative 
scheme. We illustrate the method by constructing a num- 
ber of self-consistent Kuzmin-Kutuzov models with thick 
tubes. In Section 4 we consider special and limiting cases 
for which the algorithm simplifies, and where further insight 
in the method can be gained by analytic means. Concluding 
remarks follow in Section 5. 



2 OBLATE GALAXY MODELS 

We first summarize the basic properties of oblate Stackel 
models, present the fundamental integral equation for their 
phase-space distribution functions, and then outline an it- 
erative scheme for its solution. Derivations and further in- 
formation can be found in dZ85, DZ, and ZH. 



2.1 Orbits and integrals of motion 

The motion in an oblate galaxy with a gravitational poten- 
tial of Stackel form is best described in prolate spheroidal 



coordinates (A, v, <j>). These are related to standard cylindri- 
cal coordinates (R,z,(f>) by 



R 



2 (X + a)(v + a) 



2 (A + 7 )(i/ + 7 ) 



(2.1) 



(a - 7 ) ( 7 -«) 

where a and 7 are constants and the coordinates v and A 
lie in the range — 7 < v < —a < A. Surfaces of constant 
A are prolate spheroids, while those of constant v are two- 
sheeted hyperboloids. The foci are located at z = ±-\/ 7 — a. 
Each set of (A, v, <j>) corresponds in general to two points 
(R, ±z,<f>). In these coordinates the potential V(X, v) takes 
the form: 



V 



(X + 1 )G(X)-(v + 1 )G(v) 



(2.2) 



where G(r) is an arbitrary smooth function that determines 
the shape of the potential, and r = A, v. 

The equations of motion separate in the (A, v, <j>) coordi- 
nates. Since the potential is axisymmetric, the momentum 
p<p conjugate to <f> is constant, and equals L z = R 2 (t>, the 
component of the angular momentum parallel to the 2-axis. 
The motion in A and v, i.e., in the meridional plane, is de- 
scribed by 

B(r) 



2(r + a) 2 (r + 7 ) 



(r = X,v), 



(2.3) 



fhere 



B{r) = (r + a)(r + 7 )£ - (r + 7 )/ 2 - (r + a)h - U(r). (2.4) 
and 



U(r) = _(r + a)(r + 7 )G(r). 



(2.5) 



Here E is the total orbital energy, I2 = \L 2 Z , and ^3 is the 
third isolating integral of motion given by (cf. eq. [2.13] of 
DZ) 

7 3 = hdl + Ll) + ( 7 -«) . (2 . 6) 

Each set of values of E, I2 > and ^3 for which p\ > and 
p 2 , > in some range of A and v , respectively, corresponds to 
an orbit. It is bound when E < 0. In this case the function 
B(t) generally has three zeroes for r = vo, Ai, A2, and each 
orbit fills an area in the meridional plane defined by 



- 7 < V < Vo, 



Ai < A < A 2 



(2.7) 



Orbits of this shape are usually referred to as short-axis 
tubes. 

The constants fo,Ai,A2 are functions of E, 12,1s, and 
are called the turning points of the orbit. ZH have shown 
that the relations between the standard integrals (E, 12,1s) 
and (fo,Ai,A2) can be written as 

E = U[v ,X 1 ,X 2 ], 

I 2 = ( " a " t/o)(Al+a)(A2 + a V [-«,,o,A 1 ,A 2 ], 



/ 3 



j — a 

(t/ + 7 )(Ai+ 7 )(A 2 + 7 ) 



(2.8) 



U[-y, v , Ai, A 2 ], 



j — a 

where the square brackets indicate divided differences of the 
function U(r) defined in equation (2.5). These are defined 
iteratively by 

Uir^-U^) 



U[ti,t 2 ] 



n - T 2 



(2.9a) 
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Figure 1. The three-dimensional volume filled by a short-axis 
tube orbit in an oblate Stackel model. The thin solid lines are 
the intersections of the prolate spheroidal coordinates (A, v, 0) in 
which the motion separates with the equatorial plane z = and 
with two meridional planes (R, z) at </> = and </> = tt/2. The dot 
indicates the location of the focus along the positive ,2-axis. The 
orbital volume is bounded by four prolate spheroidal coordinate 
surfaces: the top and bottom surfaces are parts of hyperboloids 
of revolution, labelled by the turning point vo, while the inner 
and outer boundaries are spheroids of revolution labelled by the 
turning points Ai and A2. 

and 

U[r 1 ,r 2 ,...,r n ]= U[Tl ' T3 '-' T " ] - U[T2 ' T3 '-' T " ] . (2.9b) 

Tl — T 2 

The ordering of the arguments is not significant. With this 
notation the function B(r) of equation (2.4) becomes 

5(r;i/ ,Ai,A 2 ) = (r-i/ )(r-Ai)(A 2 -r)[/Tr, v , \i , A 2 ].(2.10) 

Centrally concentrated models have U'"(r) > for r > — 7. 
This guarantees that all third order divided differences 
U[t\, r 2 , T3, Ti] are strictly positive (Hunter & de Zeeuw 
1992), so that B(r) has no more than three zeroes: all orbits 
are short-axis tubes. 

The three action integrals J T = (27r) -1 <j> p T dr can be 
written as follows: 



a/2 / / ■B(A;i/ ,Ai,A 2 ) dA 

A + 7 A + a ' 



J<t> = L z 



-2B(-a;v , Ai , A 2 ) 
7 — a ' 



(2.11) 



j _ _\/2 f /B(v;vo,Xi, A 2 ) dv 



v + 7 ( — a — v) 



The integrals for J\ and J v generally need to be evaluated 
numerically. 

2.2 Distribution functions 

The fundamental integral equation for the phase-space dis- 
tribution function / sm (A, v, v\, v^,, v v ) that gives rise to a 
density p m (A, v) in a gravitational potential V(X, v) is 



Pm(A, v) 



/ sm (A, v, v\, Vc/,, v v ) dvxdv^dvv. 



(2.12) 



Because V is here of Stackel form, each orbit has three ex- 
act isolating integrals of motion, so that Jeans' theorem is 
valid: / sm is a function of the three integrals of motion, so we 
can consider / sm = f srn (E, I 2 , Is), or / sm = / sm (i / o, Ai, A 2 ), 
or / sm = f S m( J\, Jcf>, Jv)- In each case, transformation of 
dv\dv<f,dv v to dEd^dls, etc., with the appropriate Jacobian 
determinant, gives the relevant form of the fundamental in- 
tegral equation (2.12). DZ (eq. [3.2]) write (2.12) in terms of 
(E, I2, 13), while ZH discuss its forms in terms of the turn- 
ing points (1/0, Ai , A 2 ) and the actions (J\, J4,, J v ) (their eqs. 
[2.23] and [2.47]). 

Since / sm is a function of three arguments, and p m de- 
pends on only two variables, many different / sm 's will be 
consistent with the same p m , so that equation (2.12) has 
many solutions. Two of these are readily available. The 
first is the special model with / sm = f sm (E,L z ), in which 
the orbits are populated such that there is no net dependence 
on the third integral. Its distribution function can be found 
by application of the Hunter & Qian (1993) method, which 
requires (numerical) evaluation of a contour integral. The 
second is the so-called thin-orbit model, in which the stars 
occupy only the short axis tubes that have no A-excursion, 
and hence lie on prolate spheroidal shells. In this case 
/sm = ftsm( Jcf>, Jv)6( J\). Bishop (1987) and ZH have shown 
that /tsm can be found by a single real quadrature. 

We are interested in distribution functions that popu- 
late not only the thin orbits, but also those with a finite 
A-extent. Instead of the turning points Ai and A 2 , we em- 
ploy the quantities 

A m = |(Ai+A 2 ), e = |(A 2 - Ai). (2.13) 

Here e > controls the 'thickness' of the short-axis tube, 
and A m indicates its mean location in the radial direction 
(Figure 1). When e = the two radial turning points Ai 
and A 2 coincide, so that the 'radial' action J\ = 0, and 
the orbit is a thin short-axis tube. The relations between 
the standard integrals (E, 12,1s) and (fo, A m ,e) follow from 
equation (2.8), upon substitution of Ai = A m — e, A 2 = 
A m + e. 

With the definitions (2.13), equation (2.23) of ZH can 
be transformed to the fundamental integral equation in 
terms of the three integrals vq, A m , e: 

p m (A, v) = 



-a co 



A m + a 



4^/d,o / dA m / de , ^(Wo.A!,^) 



yj{v -v){-a-v ){\-v ) (2.14) 



v i(A+c) |A-A m | 

c[(A m — vp) 2 — e 2 ]/ sm (t / o , A m , (.) 

A/[(A m -!/) 2 -e 2 ][(A m + a)2- e 2][ e 2_( A _ Am ) 2 ]' 



where 



U 



*_ U[vo, Ai, Ai, \ 2 ]U[vo, Ai, A 2 , \ 2 ]U[vo, vq, Ai, A 2 ] 



(2.15) 



'U[vo, Ai , A, A 2 ]?7[i/, i/o, Ai , A 2 ]?7[i/o, —a, \\ , A 2 ] 

and we still have to substitute Ai = A m — e, A 2 = A m + e. 
The area of integration in the (A m , e) plane is illustrated in 
Figure 2a. 

2.3 Iterative scheme 

Our aim is to construct distribution functions / sm (i/o, A m , e) 
that also populate orbits with non-zero thickness e > 0. 
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Figure 2. Integration areas for the fundamental integral equation, a) In the (e, A m )-plane, defined in equation (2.13). b) In the 
(s,t)-plane, defined in equations (3.2) and (3.9). The light shaded regions indicate the full integration areas: all orbits with inner 
and outer A-turning points that correspond to values of (e, A m ) or (s,t) in these areas contribute density on the spheroidal shell with 
coordinate A in configuration space. The specific choice (3.6) for the function g sm only populates orbits up to a maximum relative 
thickness s m ax, so that the integration over s runs between and s m ax, as indicated by the dark shaded regions. The thin-orbit model 
has %ax = 0, so that the integration areas shrink to a point, indicated by the filled squares. 



In the spirit of Bishop (1986), we resolve the distribution 
function / sm as the following product 



/smf^O, A m , e) — fgsm(v0, A m )^ S m(i / 0, A m , e), 



(2.16) 



where c/ sm gives the distribution of the radial excursions e of 
the orbits, which may depend on the values of the latitudinal 
turning points i>o and the mean radial positions A m . Writing 
/ sm in this way does not imply any restrictions; all distri- 
bution functions can be split up as in (2.16). If we specify 
c/sm, and substitute the result in equation (2.14), we are left 
with an integral equation for /g Sm (i / o, A m ). The choice of 
c/sm determines what function will be found. 

We choose to normalize the function c/ S m such that 



c/sm(fo, A m , e) dJx = 1, 



(2.17) 



for fixed (i/o, A m ). It then follows that when c/ S m = 8(J\), the 
function /g Sm is identical to /tsm, the thin orbit distribution 
function of Bishop (1987) and ZH. It is given by 

1 



/tsm( m, o) 87T 2 A /Am+7(Am-I / o)^[l / 0,A m ,A m ,A m ] 

x (A m + a)/9(A m ,-a)-^/(-a-fo)^[i / o,-«,A m ,A m ] 



[d(A m — <r)p(A m , a)/da]da 



(2.18) 



^/(a-i/ )U[a, f , A m , A m ] 



where p equals p m , the model density. The divided differ- 
ences of U are always derived from the model potential V, 
but the above expression gives /tsm for any axisymmetric 
density p in the potential V(\,v). 

When the function c/ S m is sharply peaked near e = 
(i.e., near J\ = 0), the solution /g Sm of (2.14) should be 
very similar to the thin orbit function /tsm- This suggests 



the following approach. We start with the zeroth-order dis- 
tribution function 



fo(vo, A m , f) = /tsm{/5m}ffsm(l / 0, A m , (:), 



(2.19) 



where /t sm {/>m} is short-hand for the thin-orbit function 
/tsm(A m , fo) that follows from equation (2.18) upon substi- 
tution of p = pm- Since c/ S m is not the delta function in J\ 
that is appropriate for /tsm{pm}, the residual density 



Pi 



/tsm{/5m}ffsm(l / 0, A m , (:) d 3 



(2.20) 



does not vanish everywhere, although we expect it to be 
much smaller than p m . We have /g Sm = fo + fc, where f c is 
the solution of 



Pi 



/c(fo, A (i/o, A m , e) d v. 



(2.21) 



This is the same integral equation as (2.14), but now for the 
density p\ in the potential V{\,v). For sharply-peaked </sm 
we approximate f c by /i = /tsm{pi}. Taking as first order 
approximation /g Sm = ftsm{pm} + /tsm{pi} then leads to a 
residual density p2 , which should be smaller than p\ . We 
can repeat this process as many times as we want, which 
leads to the following algorithm: 



/gsm(j / 0, A m ) — /tsm{p 8 }, 



where 

Po = p m (A, v), 

Pi+i = Pi 



/tsm{p 8 }</sm d 3 V, (i = 1, 



(2.22) 



(2.23) 



If the residual densities p t decrease with increasing i, the se- 
ries (2.22) will provide an increasingly better approximation 
to the actual distribution function. If this process converges, 
we still have to check that the resulting /g Sm is non-negative 
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everywhere. If it is not, it is not a physical distribution 
function, and another choice needs to be made for <7 S m- 

When c/ sm is sharply peaked, the zeroth order approx- 
imation (2.19) may already be adequate. Shu (1969) used 
it to construct self-consistent flat circular disks with nearly 
circular orbits. For less sharply peaked functions the itera- 
tive scheme (2.22)-(2.23) should work very well. However, 
as the thickness of the populated orbits increases, the den- 
sity residuals p % will become larger and it is not clear a priori 
whether the algorithm will converge. We have implemented 
the algorithm, and have constructed a number of models. It 
turns out that convergence is reached easily for models with 
quite 'fat' orbits e ~ 0.7(A m + a), but that the number of re- 
quired iterations increases strongly for broad c/ sm functions. 
The algorithm is described in detail in Section 3. 

2.4 Kuzmin-Kutuzov mass models 

We illustrate our method by applying it to the construction 
of self-consistent models with potential 

GM 



and associated density 



Pm(A, 



Mj_ a(X + 3y\v + v) 

47T 



Xv 



(XufiHVx + ^f 



(2.24) 



(2.25) 



This mass model was introduced by Kuzmin (1956) and 
connects Kuzmin's (1953) flat circular disk (7 = 0) with 
Henon's (1959) spherical isochrone (7 = a). The surfaces 
of constant density are s moo th and nearly oblate spheroidal 
with an axis ratio ~ \J 7/ 'a, and become slightly less flat- 
tened at large radii. Kuzmin & Kutuzov (1962) showed 
that the distribution function f sm (E,L z ) could be found as 
a series expansion in powers of E and L z . Many proper- 
ties of these Kuzmin-Kutuzov models were described by 
DZ, who also derived a closed form for f sm (E,L z ), al- 
beit with a typographical error (see Batsleer & Dejonghe 
1993). ZH showed that the thin-orbit distribution function 
/sm = /tsm(i / o, X m )8( J\) (Figure 3) can be given in terms of 
elementary functions, and discussed its properties in detail. 



3 DESCRIPTION OF THE METHOD 

We first discuss a practical way to choose the function c/ sm , 
introduce more convenient variables, and show that the con- 
vergence of the iterative scheme depends mostly on the mo- 
ments of </ sm . Then we show how kinematic properties of 
the models can be calculated with little extra effort, and we 
briefly describe the numerical implementation. 

3.1 Normalization of c/ sm : the function c/ sm 

The normalization of c/ sm is to some extent arbitrary, as is 
the factorization (2.16) of the distribution function. All we 
require is that c/ sm is normalized such that in the limit of 
thin orbits only, we recover /t im . The most natural nor- 
malization is to take the condition (2.17), but at constant 
values of the actions ( J$, J v ) rather than at constant val- 
ues of the integrals (fo, A m ). However, the transformation 
(vq, X m , e) <-> ( Ja, J4,, J v ) generally requires numerical inver- 
sion, so that in fact both these normalizations are not very 



practical. We therefore work with a function g sm (vo, A m , s) 
defined by 

- / -v \ c g{yo, Am) , . s /o 1 \ 

g S m(vo, Xm, e) = — - g S m(vo, A m , s), (3.1) 



m 



where 



s = - , (3.2) 

so that s is an integral of motion which gives the relative 
thickness of the orbit, and < s < 1. We require that 



#sm(fO, Xm, S) As 



(3.3) 



at fixed vo and A m . By comparison of equations (3.3) and 
(2.17) it then follows that 



Cg(vo, Xm) 



(Xm+Ot)\/X 



ds 2 tfsmf^o, A m , s)\dJ\/ds 2 \(v 0) x m ) 



(3.4) 



which therefore depends on the choice of g sm . The partial 
derivative of Ja can be written as a single quadrature, and 
is given in Appendix A. It is even in s, and hence a function 
ofs 2 . 

We have not incorporated the factor (A m -(-ci)-\/A m — pq 
in the definition of c g , because it diverges on orbits with 
X m = vo = —ce, i.e., on the oscillations along the 2-axis that 
just reach the foci of the spheroidal coordinates. We will 
come back to the behaviour of the algorithm near the foci 
in Section 4.2. In the limit g sm = 8( J\), so that only thin 
orbits are populated, we have s = 0, c/ sm (s) = <*>(s 2 ), and 



Cg^O, Xm) 



\/2(A m + 7) 

\JU\yo, X m , X m , X m ] 



(3.5) 



This is well-behaved for all physical values of — 7 < vo < 
— a < X m - 

3.2 Choice of c/ sm ; moments 

Any distribution function / sm can be written as / gim c 9 j sm . 
As an example, Figure 4 shows both the factor / gsm Cg and 
the factor c/ sm for the two-integral Kuzmin-Kutuzov model 
with f sm (E,L z ). Both factors vary smoothly, but the func- 
tion g srn shows a range of behaviour as a function of s 2 , de- 
pending on the values of A m , vq. In this paper we consider 
functions c/ sm that are even in s, and are of the form 



for < s < 



(3.6) 







with q > — 1 and < s n 



for s > s max , 

< 1. These functions are all 



normalized as in equation (3.3), and show a similar range of 
shapes as seen in Figure 4. In principle, we can choose the 
maximum relative thickness s max to be a function of vo and 
Xm, but we do not do so here, and from now on we suppress 
the dependence of c/ sm on vo and A m in the expressions that 
follow. Figure 5 illustrates the cases q = 0, 1, 2. 



We define the moments {s ™(/ S m} of c/ sm by 



1 2n \ 



2n I \ j 2 

s g srn (s)ds 



(3.7) 
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Figure 3. The distribution function /tsm(A m , vq) for an E5 Kuzmin— Kutuzov model with a. = —1 and 7 = —0.25. The thick solid 
curves are contours spaced logarithmically at intervals of 2. The focal corner lies at A m = vq = 1. See also Figure 2 of ZH. 



Figure 4. The distribution function f sm (E,L^) for an E5 Kuzmin— Kutuzov model with a = —4/9 and 7 = —1/9. The surface is the 
/g Sm part, which includes the normalizations factor c g . The thick solid curves are contours spaced logarithmically at intervals of 2. The 
focal corner lies at A m = vq = 4/9. The small diagrams in the upper half show the behaviour of g sm as a function of s, on the interval 
[0,1]. 



for n = 0, 1, . . .. It follows that (</ S m} = (s° g S m) = 1, by 
the normalization (3.3). The higher moments can be given 
explicitly for the choice (3.6): 

/ 2n i T(g + 2) 2 n /„ „\ 

( s * m >=r(n + g + 2) W (3 - 8) 



When s max — ► 0, we have c/ sm (s) = <*>(s 2 ), and all higher 
order moments vanish. When < s max <C 1, the higher 
moments decrease very rapidly with increasing n. 
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3.4 Convergence 

When s max <C 1, the function c/ sm is sharply peaked, and 
it is useful to expand it in a series of derivatives of delta 
functions (e.g., Fridman & Polyachenko 1984, p. 150): 



i« = £(-ir< 



n i 2n 



(3.12) 



where the {s 2n g srn ) cx s^ ax are the moments of c/ sm defined 
in equation (3.7), and S^™- 1 is the n th derivative of the delta 
function, which satisfies the relation 



Figure 5. Three different functions g sm of the form (3.6) with 
q = (solid line), q = 1 (dotted line), q = 2 (dashed line). They 
are normalized with respect to s 2 (see eq. [3.3]). 



3.3 New variables 

The fundamental integral equation appears in our iterative 
scheme in the form (2.23). We transform it to a more useful 
form by means of an alternative set of variables. In addition 
to the relative orbital thickness s, defined in equation (3.2), 
we introduce 



t : 



A — A n 



vo 



A 



(3.9) 



A m + a —a — v 

so that < t, < u, and < x < 1. These definitions mix 
the turning point variables (fo,A m ) with the coordinates 
(A, i/), but they facilitate the analysis of the fundamental 
integral equation. Carrying out the substitutions (3.9), and 
rearranging various terms results in 



S (n) (x)h(x)dx = (-l) n h (n) (0). 



(3.13) 



We use this expansion to show that the convergence of our 
iterative scheme depends mostly on the moments of c/ sm , and 
less on its detailed functional behaviour. 

Substitution of (3.12) in equation (3.10) gives 



Pi+i = Pi 



du wi(u) y^(s 2 "ff sm ) 

n 



2n .d n W 2 (0,u) 



d(s 2 



where we have written the t-integral as 



W 2 (s,u) 



di^M^/ tsm 



■ t 2 



{Pi}, 



(3.14) 



(3.15) 



and we have used the definition (3.13) to carry out the s- 
integration. The first term in the series expansion for c/ sm 
is <*>(s 2 ). Its contribution to the right hand side of equation 
(3.14) equals —p%. Upon substitution of the specific form 
(3.8) for the moments, we are therefore left with 



Pi+i = Pi- duwi(u) ids 2 g sm (s)x 



(3.10) 



^ w 2 (s,t,u) 



vhere we have written 
4a/2 



Wl 



W 2 



and 



'u(l — u)(l — x + xu) 
[[l-x + xu(l + t)f - 



(1- 



^/(l + txf-(l-xfs 2 ^/l- 
x U*(s,t,u)c g (t, u) 

(i+t) 3 / 2 vn? ' 



(3.11) 



■x + xu(l + t) 



and U* are defined in equations (3.4) and (2.15), 
respectively. The thin orbit function ftsm{pi} is defined in 
equation (2.18), and is independent of s. Both w\ and w 2 
depend also on the coordinates A and v, but we suppress 
them as arguments because they are not integration vari- 
ables. The square root of s 2 — t 2 vanishes at s = t = 0, 
which lies in the area of integration (Figure 2b). For this 
reason we have not incorporated it in the definition of w 2 . 
Finally we remark that w 2 is an even function of s, and 
hence depends on s 2 . 



Pi+i 



E 



r(g + 2) 

r(« + g + 2) 



du wi (u 



d n W 2 (0,u) 
d{s 2 Y 



(3.16) 



This is valid for all functions c/ sm of the form (3.6) with 

%ax ^ I- 

Since < t < s < s max , it follows that both s and t are 
small when c/ sm is sharply peaked. The functions w 2 (s,t,u) 
and ftsm{pi} then vary little over the integration area, and 
so does W 2 (s,u). Its derivatives with respect to s 2 are fi- 
nite at s = 0. Since W 2 contains the thin orbit function 
/tsm{/5 8 } as a factor, and since this is independent of s and 
proportional to p t , it follows that all the derivatives of W 2 
are similarly proportional to p t . 

For small s max , the first term in the series on the right 
hand side of equation (3.16) dominates. This means that the 
residual density p;+i is roughly proportional to the residual 
p t times the first moment of the function g sm . Since this is 
proportional to s^ ax , this shows why even for moderately 
peaked functions our iterative scheme converges rapidly. It 
furthermore shows that the shape of c/ sm is less important 
than its moments. We can vary the shape of c/ sm without 
significantly affecting the residual density, as long as we do 
not dramatically change the lowest order moments of g sm . 
Since ftsm{pi} is derived from these residuals, the solutions 
/gsm(A m , vo) that correspond to different c/ sm will be rather 
similar as long as the first moments of these c/ sm functions 
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Figure 6. Successive steps of the iterative scheme to obtain jgsm for an E5 Kuzmin— Kutuzov model with a = —1,7 = —0.25, and a 
function g sm with parameters q = 2 and s^ ax = 0.9. Shown is the residual density Pi/p m for each step. The bottom right plot shows 
the ratio of the resulting /gsm (after five iterations) and /tsm- 



Three-integral oblate galaxy models 9 



are the same. 

We found that even for very broad c/ sm functions the 
iterative scheme performed well. Figure 6 shows the den- 
sity residuals in the construction of an E5 Kuzmin-Kutuzov 
model with q = 2 and s max = 0.9. The computation was 
continued until |p;/p m | < 10 -3 , which occurred after five 
iterative steps. Even for this 'fat' model the residuals de- 
crease rapidly. The resulting / gS m is also shown, compared 
to the thin-orbit distribution function /t im . 

3.5 Velocity moments 

When a distribution function / sm = / gim c g j sm for a given 
density p m in a potential V has been found, we are most 
often also interested in the observable kinematical properties 
of the resulting dynamical model. These follow by taking 
the appropriate velocity moments of / sm , followed by a line- 
of-sight integration. It turns out that the intrinsic velocity 
moments can be found easily by a slight extension of our 
algorithm. 

The intrinsic velocity moments are defined as 



(3.17) 



The expressions for the velocity components at a point (A, v) 
along an orbit with turning points (fo,Ai,A2) are given in 
equation (3.2) of ZH. Upon transformation to our variables 
(s, t, u, x) we find 

v\ =2{\ + a) (l - X + XW \ s 2 -t 2 )U[\, A 1; A 2 , „ ], 



(1 + *) 



m(1- 



v\ =2(A + a) | U[-a, Ai, A 2 , vo], 



(1 + *) 



(3.18) 



2 .(l + xi) -(l-x)V . xrrr , , 

v 2 v =2(A-t/) 1 ^ ; \ ' — (1-«)C/[j/,A 1 ,A 2 ,j/ ]. 



(1 + t) 2 

Hence, if we insert v'^v^v* in our equation (3.10), after use 
of the transformation (3.18), we find the contribution to 
the required moment. Since the [/-functions occur in U* , 
they have already been evaluated, so calculating the velocity 
moments in parallel with carrying out the iteration to get 
/ sm adds very little to the required CPU time. 

As an example, we have constructed two E5 Kuzmin- 
Kutuzov models with q = and s max = 0.5 and 0.7, re- 
spectively. In five iterations the residual density is less than 
10 -3 /9 m . The velocity moments are shown in Figure 7 (cf. 
Figures 5 and 7 in ZH). As expected, the dispersion in the 
v- and (^-directions decrease while the 'radial' dispersion 
increases. The s max = 0.7 model has the largest ratio of 

/ '( v 4>) ) which is of the order of 0.25. 



3.6 Numerical implementation 

We have written a code to implement the iterative scheme 
defined in equations (2.22) and (2.23). The density residu- 
als and the terms in the series for / gS m are calculated on a 
(A, v) grid that doubles as a (A m ,i/o) grid. The quantities 
are interpolated using splines as their values are also needed 
in between mesh points to evaluate the s- and t-integrals. 
The actual integrations are carried out using Romberg inte- 
gration after switching to more appropriate integration vari- 
ables instead of s and u to remove the square roots from the 



denominators. 

The (A, v) grid points have to be chosen with some care. 
We need to cover the full A domain in order to prevent 
boundary errors in the determination of the residual den- 
sity. This is accomplished by using a grid that is linear in 
the variable t\ defined by 



X(tx) + a = (Ajj + a) 



(2-2t x )p" 



< t x < 1, 



(3.19) 



where An is the value of A in the middle of the grid, p 
determines the resolution near A = — a and p is chosen to 
match the large-radii behaviour of the distribution function. 
Most of the computations in this paper were done using 
p = 3, p = 0.25 and \ H = 4. 

When the model is very flattened towards the equatorial 
plane, its associated prolate spheroidal coordinate system is 
very elongated along the 2-axis. The net effect of this oppo- 
site elongation is that the density and distribution function 
are sharply peaked near v = —j. We therefore use a i/-grid 
that is linear in the variable u v , where 



f(« I /) + 7 = A 



1 



(3.20) 



.B + il + u^P- 1 5 + 1. 

for some p > 1; B is set to (\) v ~ 1 ■ The parameter A is 
determined so that J'(l) = —a. We have found this substi- 
tution to be adequate for models as flattened as E7. 

The thin-orbit distribution function is computed from 
the density residuals using (2.18). A straightforward imple- 
mentation of (3.10) is feasible, but a single iteration step on 
a (A, v) grid of 50x50 would take about 100 minutes CPU 
time on an HP735. There is a faster way to implement 
the algorithm. The only part of the integral in (3.10) that 
changes in successive iterations is the /tsm{/5 8 } factor, which 
does not depend on s. We therefore exchange the order of 
integration, and carry out the s-integral first. It is 



Tg(t, U, X) 



ds 



2 w 2 (s,t, u)g srn (s) 
x/s 2 - t 2 



(3.21) 



and can be tabulated before the first iteration. Using a 
81x27x100 grid for (t,u,x), the initialization stage takes 5 
minutes CPU time, or 10 minutes if the velocity moments 
have to be computed as well. Each iteration step then sim- 
plifies to 



Pi+i = Pi - du wi(u) / dt T g (t,u, x)f tsm {p t } 



(3.22) 



which takes about 1 CPU minute to complete. The program 
is written in Fortran. The number of iteration steps to reach 
a relative accuracy better than 10 -3 is 1 for s max = 0.1, 3 
for s max = 0.5 and 5 for s max = 0.7 in the case of q = 
models. It is smaller when q > 0. 

There is a further reduction possible in the limiting 
cases of a spherical or a disk galaxy. These are described 
in Sections 4.4 and 4.5, respectively. 



4 SPECIAL CASES 

Approximate solutions of the fundamental integral equation 
(3.10) can be found by analytic means at large radii, and 
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Figure 7. The intrinsic velocity moments for the E5 Kuzmin— Kutuzov thin-orbit model (solid) and q = models with s max = 0.5 
(dotted) and 0.7 (dashed lines). 



near the foci of the spheroidal coordinates. We consider 
them in turn, and then show how the scheme simplifies for 
models with sharply peaked c/ sm functions, and in the spher- 
ical and disk limits. 



4.1 Large radii 

The relation (3.10) for the density residuals simplifies con- 
siderably in the limit of large radii, i.e., A — > oo so that 
x — > 0. The functions w\ and W2 then reduce to 



wi (u) 



W2(s, t, u) 



4a/2 



\A(i-**)' 

U*(s, t, u)c g (t, u) 

(1 + t) 3 / 2 • 



(4.1) 



In order to evaluate the various third-order divided differ- 
ences that appear in the definitions (2.15) and (3.4) for U* 
and c g , respectively, we assume that s max < 1, so that 
the orbits that contribute to the density at (A, v) have 
A 2 > A > Ai = A m (l - s) > A m (l - s max ) > -a. The 
details of the calculations are given in Appendix B. We find 
that in this case neither U* nor c g depends on u at large A, 
and that 



U*c 



GM (1 + t) 
2 5 / 2 C A 



L(s,t), 



(4.2) 



where L(s,t) is the elementary function given in equation 
(B4), and the constant C g > 1 is defined in equation (B10). 



C a 



1 in the thin orbit limit g sr , 



S(s 2 



We consider a flattened density p, that falls off as a 
power of A at large radii: 



ft(A, v) 



(4.3) 



for some (positive) value of p. Selfconsistent models with 
finite total mass must have p > 3/2. By Kuzmin's formula, 
such models must also have p < 2 (e.g., de Zeeuw, Peletier 
& Franx 1986). Non-consistent densities p m may fall off 
steeper than this. Substitution of the form (4.3) in expres- 
sion (2.18) for /tsm{ft} and use of the approximations (Bl) 
then shows that /tsm oc \\^ p times a function of vo (cf ZH, 



eq. [2.56]). Transformation to the variables (s,t,u) gives 



/tsm{p 8 } 

where 



Xp- 1 ir 2 GM 



fura = ft (-«) + a/« 



[dp^(«i')/(2«i'] du' 
t/m — u' 



(4.4) 



(4.5) 



Thus, at large radii the thin orbit distribution function be- 
comes a product of a function of t and a function of u. 

Substitution of all the above approximations in the ba- 
sic relation (3.10) shows that for A ^> —a the triple inte- 
gration over s,t and u reduces to the product of an integral 
over u times a double integral over s and t: 



p,+i ~ p, 



Lg(p) 1 f du ft sm (u) 



*Cg * P J ^/u(l-u) 



where we have defined 



L g(p) = - ds Ssm(s) 

7T 



dtL(s,t) 



(1 +t)2-P^fs 2 ~^¥ 



(4.6) 



(4.7) 



The M-integration in equation (4.6) can be carried out, and 
equals ■KpX(v). Since L g (jp) is a constant, it follows that the 
triple integration over /tsm {ft} is proportional to p^(p)/\ v , 
i.e., it is proportional to p t itself: 

LAP)] 



p t+ i ~ p % 1 



C 



(4.8) 



Thus, the residual density at large radii becomes smaller by 
a constant factor [1 — L g (p)/C g ] at each iteration step. 

In the thin orbit limit c/ sm = <*>(s 2 ), and we have 
C g = Lg(jp) = 1, so that Pi-yi = for all i. This is as it 
should be, since / gS m equals /t S m{ft-n} exactly in this case. 
Equation (4.8) implies that for broadened functions c/ sm the 
distribution function / gS m at large values of A m is given by 

C 

/gsm(>0,A m )~ 9 /tsm{Pm}, (A m >-£*), (4.9) 

L g(P) 

so that we can find it without iteration. 
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Figure 8. Limiting behaviour at large radii of various properties of models with a density fall-off p m oc l/A 2 , and functions g sm with 
q = (solid), 1 (dotted), and 2 (dashed), and values of s m ax between and 1. Here A ~ r 2 , and /gsm becomes a constant factor 
C g /Lg(2) times the thin-orbit distribution function /tsm (see eq. [4.9]). Shown are, as a function of s m ax: a) the factor C g /L g (2), b) 
the intrinsic mean streaming motion {v^) in units of the maximum possible streaming (w^)tsm, assuming all stars have the same sense of 
rotation around the symmetry axis, c) the second moment (v 2 ^ ) of the 'radial' velocity, in units of GM/s/X, and d) the second moments 
(vV) and (v 2 ) in units of their values in the thin-orbit model. These results are valid for A m (l — s m ax) ^ —a. 



The values of the constants C g and L g {jp) can be found 
by numerical evaluation of the integrals (BIO) and (4.7). For 
sharply peaked g sm they can be approximated by expanding 
the integrands in powers of s 2 , and evaluating term by term. 
This gives 



, / 1005 _£71 ,227 2_n 3 J_ 4w 4 v 
+ 1 128 128-^+ 128-P 32^+32^ A* # sm / 

+ 0({s 6 g sra )), 

so that 

Pi + l - Pi [(^-|-P+ip 2 )( s2 5 , sm> 

+(i-ffp+if/-i/+^ 4 )(^ m >+...]. 



(4.10) 



(4.11) 



This shows again that for sharply peaked c/ sm its first mo- 
ment is mostly responsible for the convergence of the itera- 
tive scheme. 

In a similar fashion it is possible to derive approxima- 
tions for the intrinsic velocity moments. Substituting ap- 



proximations (Bl) in (3.18), we find that the velocities can 
be written as a 

, GM 



e T (u)Ll(s,t), (r = \,v,<t>), 



(4.12) 



where L T V is given in (B12), and l T = 1,(1 — u),u for 
t = \,v,<f>, respectively. Substitution of (4.12) and (4.9) 
in (3.17) yields 



Pm 



,^ n) _ L™ 1 /GMW 2 I Au fl m {u)i T {u) (4i3) 



L g (p) tt\p V VX / J J 





phere we have introduced 



L T g "(p) = - ds 2 g sm (s) 



AtL T v {s,t)2 L{s,t) 



(4.14) 



Again, the M-integration can be carried out. Since £J(0, 0) = 
1 for t = v,4>, the M-integral equals p m (f")tsm and 
Pm("J}tsm, respectively, which are the velocity moments of 
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the thin-orbit model. For r = A the M-integral equals p n 
hence 



L g ( P ) ' 



(r = v, <t>), 



(4.15) 



For sharply peaked c/ sm functions the dispersions and rota- 
tion velocity can therefore be approximated by 

W) ^ (»*)t S m [l + (If - \p){s 2 g,m) 

+ (irk -mp + ikp 2 )^^) + ..]; 
M>^^f[l< s 2 , sm > + (^-lp)< s V m > + ...]; 
<^> ~ <^> tsm [1 + (-f + f p)(s 2 g sm ) (4.16) 

■ / 531 215 3 2w 4 \ . 1 

+ (-128 ~ — p ~ sP )( s J-) + •••]; 

(«|) ~ (vfytsm [l + (| - ip)(s 2 filsm> 

+ (i%-||P+|/)^ 4 ^ m > + ...]. 

Figure 8 illustrates the behaviour of C g /L g (p) for our choice 
of </sm, for p = 2. It also shows intrinsic velocity dispersions 
and rotational velocity in terms of the thin-orbit values. 
The above results show that for a given density distribution 
Pm(X, v), and with c/ sm a function of s only, the kinematic 
properties at large radii do depend on A and v, but they 
follow from those in the thin orbit model by multiplication 
with a factor wich depends on c/ sm and p only. 

Equation (4.8) shows that our iterative scheme will con- 
verge for any (/ S m for which < L g (p)/C g < 2. A useful 
upper limit for L g (p)/C g can be obtained by evaluating 



max 

0<i<i m „ 7Tft(s) 



dtL(s,t) 



(l + i)2 



(4.17) 



-t 2 



It is easily verified by numerical integration of (4.17) that 
L g /C g < 2 as long as p < 9/2. This includes all physically 
relevant cases, so our iterative scheme will always converge 
at large radii for all c/ sm functions. 

4.2 Behaviour near the foci 

We now investigate the behaviour of the iterative scheme 
near the foci of the spheroidal coordinates, where A = v = 
— a. The variables s,t,u and x that appear in the basic 
relation (3.10) for the density residuals take their full range 
of values near the foci, but the triple integration over s,t 
and u can nevertheless be simplified, because the factors 
/tsm{p 8 }(t, u), U*(s,t,u), and c g (t,u) which appear in the 
integrand, all simplify. 

The only orbits that contribute density at the foci in 
the thin-orbit model are those with A m = vo = —a. Or- 
bits with A m > —a lie on spheroidal shells which inter- 
sect the 2-axis above the foci, while orbits with A m = —a 
and vo < — ce are 2-axis oscillations that do not reach the 
foci. It follows that p m ( — a, —a) is determined exclusively 
by /tsm{pm}(- ce, — a). The relation is (see ZH, eq. [2.44]) 

Pm( — a, — a) [1 + x ] 



/tsm{pm}(-a,-a) 
where 



x 



-a — vo 



8ir 2 \Jj — ceU[—a, — a, — a, — a] 
xu(l + t) 



A n 



-v 



l — x-\-xu(l-\-t) 



(4.18) 



(4.19) 



and U[—a, — a, — a, — a] = U"'( — a) > 0. In the limit where 
A m = vo = —a, the value of xo can still lie anywhere between 
and 1, so that the factor in square brackets in equation 
(4.18) can take any value between 1 and 2, depending on the 
direction along which the focal corner in the (A m , i/o)-plane 
is approached. ZH refer to this property of /tsm by saying 
that it has radial behaviour in the focal corner (Figure 3, 
and ZH Figure 2). Without this behaviour of /tsm it would 
not reproduce the correct density p m ( — a, — a). 

Since the thin-orbit distribution function ftsm{pi} is 
proportional to p m ( — a, —a), the residual density p\ depends 
only on the local behaviour of the distribution function near 
the focal corner. This means that we can approximate U* in 
equation (3.10) by the constant value [U"'( — a)] 3 ^ 2 , so that 
it can be taken outside the triple integration. 

We show in Appendix C that near the focal corner the 
function c g can be approximated as 



Cg(vo, \m)- 



0(7-") 1 
y/U"'(-a) Jg(xo)' 



(4.20) 



where J g (xo) is a weighted integral of the function </ S m, de- 
fined in equation (C7). As a result, c g also has radial be- 
haviour near the focal corner, except in the thin-orbit limit 
when J g (xo) = 1. 

Upon substitution of these approximations in relation 
(3.10), we find 



pi(-a,-a) ~ p m ( — a, — a)[l - F g (x)], 



(4.21) 



where 



F g (x) 



du 

71-2 J \/u(l — u)(l — x-\-xu) 



ds 2 I At g sm (s) (1 + xq) (4.22) 

VT^J (1+t) 3 ' 2 Jg(xo) 
-s 

[[l-x + xu(l + t)] 2 -(l-x) 2 s 2 ] 
^J(l + tx) 2 -(l-x) 2 s 2 ^Jl-x + xu(l + t) 

and we still have to substitute relation (4.19) for xo = 
xo(x,u,t). When c/ sm (s) = 8(s 2 ) the integration over t and 
s gives 7r since then J g (xo) = 1, and so does the remaining 
integral over u, so that then F g (x) = 1 and the residual 
pi = 0, as it should be for the thin-orbit model. However, 
when g sm is not infinitely sharply peaked, the value of p\ at 
the foci (A = v = —a) depends on the direction of approach, 
i.e., on the value of x. Thus, the residual density has radial 
behaviour near the foci. 

The triple integral (4.22) requires numerical evaluation 
for < x < 1 and general c/ S m- We have found it to be a 
slowly varying monotonic function of x for our choice (3.6) of 
functions g sm (s) (see Figure 9). It is bounded by the values 
F g (0) and F g (l). We show in Appendix D that the triple 
integration reduces to a single integral for x = and x = 1, 
and that the remaining integrals can be found explicitly for 
our set of functions (3.6). By combining equations (C9) and 
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Figure 9. The behaviour of (1 + xo) / F g (xo) for (a) q = 2, (b) q = 1 and (c) q = 0, as calculated with our iterative scheme, for 
s max = 0, . . . , 0.9 in steps of 0.1. The factor l/F g (xo) decreases monotonically with s m ax- At most 10 iterations were computed for each 
model, regardless of the achieved accuracy. The scheme did not converge for the (q = 0,s^ ax = 0.9) model. In (d) and (e) the results 
(crosses) from (a)— (c) are compared to the analytical value for q = (solid), 1 (dotted) and 2 (dashed) as a function of s m ax- In (f) the 
lines Fg(l) = 2 (dotted) and F g (0) = 2 (solid), which are in effect a relation between q and s ma x, are plotted for our g sm functions. In 
the shaded area both F g (0) < 2 and F g (l) < 2, which is an indicator for convergence of the iterative scheme. 
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2Ji(f,f;2 + g ;3 



2 ) 

max ) 



(D6) we find, for q > 

2^(1,1; 2 + g; S Lax)' 
while for g = we have (cf. eqs [CIO] and [D8]) 
3[K(k)-(l + s m ^)E(k)] 



2 F 1 (l,l;2 + q;s 



2 ) 

max / 



(4.23) 



F g (l) 



(1 + s max ) [£(*) - (1 - s max ) K(k)] 
(1 + a/T^ 



(4.24) 



2sl 



ln(l-s max ), 



where K and E are the complete elliptic integrals of the 
first and second kinds, and k 2 = 2s max /(l + s max ). For each 
q > these functions increase monotonically with increasing 
s max - In the limit s max — > 1 each function reaches a finite 
value, 



F a (0) 



q + 



F g (i) 



g + f 



(4.25) 



provided q > 0. However, F g (0) and F g (l) diverge logarith- 
mically when q = and s max — > 1, although their ratio ap- 
proaches 3/2, in agreement with the result F g (0)/ F g (l) — > 
(g + f)/(g + which follows from equation (4.25). For 
given q and s max the total relative variation of F g (x) be- 
tween x = and x = 1 is therefore never larger than 1.5. 
Figure 9 illustrates the behaviour of F g (0) and F g (l) as a 
function of s max for q = 0, 1 and 2. 

The radial behaviour of the residual density pi at the 
foci of the spheroidal coordinates, and the fact that the dis- 
tribution function / sm = / Esm j im is determined by the local 
density distribution when s max < 1, means that / gsm (j'o, A m ) 
must have radial behaviour near the focal corner vo = X m = 
— a, so that it must be of the form 

/gsm{Pm}(-a,-a) = ftsra{Pm}{-a,-a)K g (x ) 

= p m (-a,-a)(l + x )K g (x ) (4.26) 
8tt 2 a / ; P^ U"'(-a) 

with K g a function of the variable xo defined in equation 
(4.19). Since x = corresponds to xo = 0, and since simi- 
larly x = 1 corresponds to xo = 1, it follows that 



KM 



1 



F g (0) 



K g (l) 



1 



F g (iy 



(4.27) 



This suggests — but does not prove — that K g (xo) has a 
modest variation with xo- It can be computed as follows, at 
least in principle. We insert K g (xo) as a factor in the triple 
integral on the right-hand side of equation (4.22), and put 
the left-hand side equal to 1 for < x < 1. Transformation 
of the variables (u,t) to (xo,t) then allows one to carry out 
the t and s-integrals. This leaves a one-dimensional integral 
equation for the function K g (xo)- In practice this must be 
solved numerically, and it is in fact easier to simply use our 
iterative scheme. We have applied it to compute K g (xo) for 
our functions c/ sm with q = 0, 1 and s max = 0,0.3,0.5,0.7 
and 0.9. The resulting functions are shown in Figure 9a,b,c. 
They are very nearly linear up to values of s max as large as 
0.7. This means that to good approximation we can take 



T- i \ 1 — zo , 
K g (xo) ~ „ + 



Xo 



F g (0) F g (l)' 



(4.28 



so that we can find the local behaviour of the distribution 
function near the focal corner without iteration. 

The results presented in Figure 9d,e also show that the 
jump in the value of / gS m at the focal corner is a function of 
s max - The magnitude of the jump going from xo = to xo = 
1 is 2Kg(l)/K g (0) = 2Fg(0)/F g (l), and hence varies from 2 
when s max = in the thin-orbit model to (2g + §)/(? + §■) 
when s max — > 1. 

The convergence of the iterative scheme near the foci 
of the model cannot be derived by investigation of only the 
residual density pi, as we have done in the above. However, 
equation (4.21) suggests that in cases where < F g (x) < 2 
also the higher order residuals \p t \ with i > 1 will decrease 
in size, so that the scheme very likely will converge. This 
condition on F g (x) is met for all our functions c/ sm when 
q > 3/4, and is also met for a large range of s max when 
-1 < q < 3/4 (see Figure 9f). 

The above results also hold when s max = 1 as long as 
q > 0. In this case there are orbits with arbitrarily large 
outer turning point A2 but low angular momentum that still 
provide density at the foci, but their contribution is van- 
ishingly small. However, when q = and s max = 1 this is 
no longer so, and our derivation of the approximate relation 
(4.21) is invalid. The density close to the foci depends on 
the details of the distribution of non-local orbits. 

We remark that all self-consistent oblate models with 
density p m in a potential V of the form (2.2) have distri- 
bution functions with radial behaviour in the focal corner, 
unless q < and s max = 0. From our analysis it is not 
clear whether the radial behaviour is still present in the lat- 
ter case. However, the two-integral distribution function 
fsm(E, L z ) is generally well-behaved at the focal corner (Fig- 
ure 4). When written in the form (2.16), it leads to a func- 
tion c/ sm which does not depend on s alone. At the focal 
corner it is identical to our q = 0, s max = 1 function, and 
along vo = —a it has q < behaviour. Our iterative scheme 
with /tsm as initial guess for / gS m fails to converge in this 
case. 



4.3 Models with peaked c/ sm functions 

We have seen that the convergence of the iterative scheme is 
determined by the moments of the c/ sm function. For models 
with sharply peaked c/ sm functions only the first moment is 
significant; higher moments can be neglected. In this case 
only a single iteration is required to approximate / gS m to 
high accuracy. 

The expression for the residual density p\ simplifies con- 
siderably. Since \t\ < s < s max and s max is small, the 
W2(s,t,u) function in (3.10) can be expanded in a Taylor 
series in t and s 2 : 



/ \ dw2 2 dw2 1 d W2 2 / \ 

W2(s, t, u)~W2 + — ^ 2 + -£t + f^/i 2 + . . . , (4.29) 

where the derivatives are evaluated at s = t = 0. Similarly, 
/tsm(i, u) can be written as a Taylor series in t. The s- and 
t-integrations can be carried out. The t«2(0, 0, M)/t sm (0, 0) 
term yields the model density po, hence the residual density 
Pi is 



^ Pi = {s 2 9sm)&p 



(4.30) 
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Figure 10. The function Sf/ft sm for an E5 Kuzmin— Kutuzov model as obtained from a model with s ma x = 0.1 and (a) q = 0, (b) q = 1 
and (c),(d) q = 2 (shaded surfaces). The difference of Sf/ft sm obtained from (a)— (c) s m ax = VO.l and s m ax = 0.1 and (d) q = and 
q = 2 models are shown as wireframe surfaces. 



with 



4.4 Spherical limit 



Sp-- 



du w\ (u) 



a - 9 + Adt 2 F u ™ +W2 ±dt 2 



dw 2 
Js 2 



,(4.31) 



with %»2 and /tsm evaluated at s = t = 0. Hence 5/9 is 
independent of the shape of </ sm . This relation even holds if 
c/sm also depends on A m and i>o- The resulting / gS m can be 
approximated by 



/gsm ~ /tsm + {s 2 g sm )Sf 



(4.32) 



where 5/ = ftsm{8p} is also independent of g sm . Thus, all 
models with sharply peaked g sm functions can effectively be 
described by the one-parameter family (4.32). 

We have computed 8f for an E5 Kuzmin-Kutuzov 
model with q = 0,1,2 and s max = 0.01 and 0.1. The re- 
sults are displayed in Figure 10 as shaded surfaces. To ver- 
ify that 8f is indeed independent of </ S m, i.e., of q and s max , 
we compare the difference of two <5 /-functions obtained for 
a different set of (q, s max ) with 8f itself. It is clear from 
Figure 10 that the differences are indeed much smaller, so 
that (4.32) holds to high accuracy. 



When 7 — > a, the prolate spheroidal coordinates (A, v, <f>) 
reduce to spherical coordinates (r,d,cf>) with r 2 = A + a 
(e.g., dZ). The potential (2.2) now is spherical, and equals 
V(r) = — G(A). The thin-orbit model reduces to the spher- 
ical model built exclusively with circular orbits (ZH, §IIg). 
Our iterative scheme will produce spherical models with a 
pre-assigned distribution of relative orbital thickness. We 
summarize briefly how the algorithm simplifies in this case. 

In the limit y = a, we must have v = vo = — a, so 
that x = 0. This reduces the expressions for w\ and W2 that 
appear in the basic relation (3.10) to the forms already given 
in equation (4.1). Furthermore, the thin-orbit function /tsm 
becomes a function of A m only. It can be evaluated without 
the need for integration since the density at radius r m , say, 
is provided only by the circular orbits with radius r m . We 
write /tsm as (cf. ZH, corrected for a typographical error of 
a factor of 3) 



/tsm{p 8 }(r m ) 



Pi{Tn 



It 2 Tm (^rn ) 



(4.33) 



where Kq is the epicyclic frequency, given by r«o(r) 
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3V'(r) + rV" '(r), and r m is denned by r 2 m = A m + a. Also, 

(4.34) 



S g (r m ) 

where S g is the integral in the denominator of the definition 
(3.4). It contains the factor \dJ\/ds 2 \ which is given in 
equation (A3), and must be evaluated with —j = vo = —a. 
Finally, the function U* defined in equation (2.15) reduces 
to 

U[— a, Ai , Ai , A2]?7[— cx, Ai , A 2 , A2] 



U* 



(4.35) 



^U[-a,\ 1 ,\,\ 2 ] 

We write r\ = Ai + a and r\ = A2 + cx, so that n and T2 are 
the inner and outer radius reached by an orbit, and r 2 m = 
+ rl). Since (7(A) = -(A + a) 2 G(A) = rVjr), it follows 
from the definition (2.9) of the divided differences that the 
function U* is a combination oiV(r) and its derivative V'(r) 
evaluated at r, n and r2. These are related to the variables 
s and t by 



2 2 
r 2 -ri 



t 



(4.36) 



so that U* = U*(r,s,t), c g = c g (r,t) and /tsm = /tsm^, t). 
This means that we can carry out the M-integration in equa- 
tion (3.10) to give 7r. 

Substitution of the above results and exchange of the 
order of the s- and t-integrations, then reduces the basic 
relation (3.10) for the residuals p t to 



Pi+i{r) = Pi{r) 

7T 



dt 



H g (r,t) r 2 m p l {r m ) 
(1 + i) 3 / 2 S g (r m )K 2 (r m )' 



(4.37) 



vhere r m = r(l + t) 1//2 , and we have defined 



H g (r,t), 



ds 2 



-J7 2 



-t- 



■ g S m(s)U*(r, s, t). 



(4.38) 



Hence, we can choose a function c/ sm , evaluate H g (r, t) once, 
and then compute p;+i at each radius r as a single weighted 
integral over p;(r). In practice it is convenient to use r m 
rather than t as the integration variable in equation (4.37). 
No further work is needed to find the entire distribution 
function / sm = / Esm j im of the model, because it is given by 



/sm(^T7 



T m ) f^n f T m ) ^ 



7T^ T m Sg (fm)'i'o('' , "i) 



(4.39) 



Calculation of the distribution function is thus considerably 
faster than in the flattened case, where computation of each 
p % requires a double integration over a distribution function 
/tsm which itself is evaluated as a quadrature. 

We remark that the arguments r m and s of the above 
spherical distribution function each depend on the classical 
integrals of motion E and L 2 . The relations follows from 



E : 



f 2 ,2 2 ^2)-^) 
1j — Z,T 1 T 2 2 2 • 



(4.40) 



with t\ = r 2 n (l — s) and r\ = r 2 n (l + s). For given r m and 
s these relations must generally be inverted numerically to 
give the associated E and L 2 . 



4.5 Disks 

The iterative scheme (2.22)-(2.23) can also be applied in the 
limit where the density flattens to a circular disk with sur- 
face density E m (A) = J p m dz, where A + a = R 2 . The only 
orbits that can now be populated are those in the equatorial 
plane z = 0, so that we must have v = vo = — 7 and hence 
u = 1 in our fundamental relation (3.10). This leads to a 
number of simplifications. 

When 11 I 1, the thin-orbit function /tsm, defined in 
equation (2.19) can be approximated as 



/tsm{p*} 



1 



1 



87T 2 -\/Am+7 U[-~j,\ m ,^m,^m\ 



A /y[-7,-«,A m ,A m ] ]jm j- d } Pi (\ m ,u')Au' \ (4 ' 41) 
y/U[—i,—i,\ m ,\ m ] 'nXduJ y/u-u' i' 



where p l (\ m , v.') = £;(A m )6(«i' — 1). The function c g now is 

(4.42) 



(X m + a)y/X 



with D g the integral in the denominator of equation (3.4), 
evaluated at v = vo = — 7- The function U* reduces slightly, 
and becomes 

U* = «7[-7,Ai,Ai,A 2 ]?7[-7,Ai,A2,A2] 



U[~7, -7, Ai, A 2 ] 



(4.43) 



«7[- 7 ,Ai,A,A 2 ]?7[-7,-a,Ai,A2]' 



If we write R 2 = Ai + a, R\ = A2 + ex, R 2 m = A m + a, and 
use U{\) = -(A + a)(A + 7)G(A) = R 2 (R 2 + 7 - a)V(R), 
with V(R) the potential in the equatorial plane, then we 
can express all the above third-order divided differences in 
terms of V and its derivatives. In particular, we have 

U[-y,-a,X,X] = ^Ql(R), 



2 L 

U[-y,\,\,\] = ±K 2 (R), 



(4.44) 



with Ko the epicyclic frequency defined in Section 4.3, and 
fio the circular frequency, given by i?fio(i?) = V'(R). The 
cylindrical radius R coincides with the spherical radius r in 
the equatorial plane, so equation (4.36) can be used to find 
U* = U*(R, s, t), c g = c g {R, t) and / t8m = f tam {R, t). 

We substitute the above approximations in equation 
(3.10), and integrate it over z in order to obtain the ba- 
sic relation for the residuals in the surface density. The 
M-integration can be carried out, and we are left with 

E i+ i(fl) = Ei(fl) 



dt H g (R, t)R 2 rn ^lo(R m )Y lt (R m ) 
(l + ty/^l + xty^DgiR^nKRm)' 



(4.45) 



where i? m = + 1//2 , x = (7 — a)/{R 2 m + 7 — a), and we 
have defined 



H 



ds 2 g sra (s)U*(R,s,t) 



^/(l + xt) 2 -(l-x) 2 s 2 



(4.46) 



Just as in the spherical limit, the iterative scheme (2.22)- 
(2.23) simplifies considerably. We can choose a function 
</sm(s), integrate it to get H g (R,t), and then evaluate £ 8 +i 
at radius R by a single quadrature. We have written c/ S m 
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here as a function of s alone, but a dependence on R m can 
be included easily. 

The three-dimensional distribution function of an in- 
finitesimally thin disk can be written as 



fsm(J\, J<f>, Jv) = fdisk(J\, J <p) 



2?r 



(4.47) 



where the division by 2ir ensures that fdisk(J\, J#) is the 
proper distribution function for the disk considered as a two- 
dimensional system. It follows from equations (2.22) and 
(4.45) that our scheme gives /disk as 



/disk(-Rm, s) 



flo(R m ) ffsm(s) 
WK 2 {Rm) D g {R m ) 



(4.48) 



The distribution function (4.48) can be written as a function 
of E and J$ = L z by use of equation (4.40), with r replaced 
by R, and L 2 by L\. 

We conclude that our iterative scheme provides a swift 
way to construct distribution functions for spheres and disks 
with a chosen distribution of the relative weights of orbits 
with different 'thickness' r\ — r\ and mean radial extent r 2 m . 
Based on our results in Sections 4.1 and 4.2, we expect the 
scheme to converge quickly, except for choices of c/ sm that 
put a lot of weight on radial orbits, i.e., that have c/ sm (l) > 0. 



5 CONCLUDING REMARKS 

We have presented a simple numerical scheme for the con- 
struction of three-integral distribution functions for self- 
consistent and non-consistent oblate galaxy models with a 
potential of Stackel form. The intrinsic velocity moments 
can be computed simultaneously. The algorithm allows one 
to choose in advance the distribution of the inner and outer 
turning points of the short-axis tube orbits that are popu- 
lated. It then derives the entire distribution function from 
the density distribution by means of an iterative process that 
starts from the explicitly known distribution function of the 
thin-orbit (maximum streaming) model, in which only the 
tubes with equal inner and outer turning points are occu- 
pied. We have shown that this scheme works well, and is 
capable of producing tangentially anisotropic models with a 
substantial radial velocity dispersion within a few iteration 
steps. The algorithm simplifies considerably in the spherical 
and disk limits. 

Dehnen & Gerhard (1993) have shown that three- 
integral flattened models display a large variety of observable 
kinematic properties, which include the line-of-sight mean 
velocity and velocity dispersion, as well as the entire dis- 
tribution of the line-of-sight velocity (the velocity profile), 
all as a function of projected position on the sky. The ob- 
servable kinematics of the tangentially anisotropic models 
constructed here can be computed in a straightforward way 
by numerical integration of the velocity moments and the 
distribution function, all of which are given with high accu- 
racy by the algorithm. 

We have investigated three special cases where three- 
integral distribution functions can be found without itera- 
tion. 

First, models that have modest radial dispersions can 
be approximated adequately by a one-parameter family of 
distribution functions, which is insensitive to the detailed 



shape of the assigned function c/ sm , but depends only on its 
first moment. We will use this family in a subsequent paper 
to investigate the stability of cold oblate models. 

Second, the structure of the model near the foci of the 
prolate spheroidal coordinate system in which the equations 
of motion separate provides information on the convergence 
of the algorithm. When the function c/ sm is chosen such that 
only a vanishingly small number of orbits with L z = and a 
large outer turning point are occupied, the density near the 
foci is determined locally, i.e., by stars on orbits that are very 
close to the 2-axis oscillations that just reach the foci. We 
have derived the local behaviour of the distribution function 
in all such models, and we have shown by analysis of the first 
residual density that the algorithm is very likely to converge 
in these cases, as indeed found numerically. However, when 
L z = orbits with large outer turning points contribute 
significantly to the density at the foci — which occurs in the 
f(E, Lz) model, and in strongly radially anisotropic models 
— our algorithm appears to have problems, at least when 
we take /tsm as initial guess for /g Sm - In view of Bishop's 
(1986) work, we expect that a similar iterative scheme can 
be used for such models, but with f(E, L z ) as zeroth order 
distribution function. 

Third, the distribution functions of the models also sim- 
plify at large radii. There they reduce to a known fac- 
tor times the distribution function of the thin-orbit model, 
which can be calculated easily. The internal velocity mo- 
ments similarly simplify at large radii. This is useful, as it 
allows a straightforward calculation of the observable kine- 
matic properties in the outer regions of these anisotropic 
flattened models. We intend to do so in a future paper. 
Absorption line kinematic measurements of elliptical galax- 
ies now reach beyond two effective radii, and a comparison 
of these data with anisotropic models of the kind produced 
by our algorithm should provide further constraints on the 
presence and shape of a massive dark halo and the dynamics 
of the outer luminous regions of these systems (e.g., Carollo 
et al. 1995). 

Finally, we remark that the iterative scheme is not re- 
stricted to oblate galaxy models. Prolate Stackel models 
have two families of tube orbits, and the thin-orbit solutions 
have been described by Hunter et al. (1990). By apply- 
ing our algorithm separately to the two tube orbit families, 
we can construct models with thick tubes. Triaxial Stackel 
models contain three families of tube orbits as well as box 
orbits. The thin-orbit distribution functions for all three 
tube families can be found by simple quadratures (Hunter 
& de Zeeuw 1992; Arnold, de Zeeuw & Hunter 1994), and 
these can again be thickened by our algorithm. The tube 
orbits account for part of the density; the remainder must 
be reproduced by the box orbits. Their distribution func- 
tion can be found by (numerically) solving a set of linear 
equations after the tube orbits have been populated. This 
last construction step is the same in thin and thick orbit 
models. Work on these triaxial models is in progress. 



It is a pleasure to thank Richard Arnold and Ma- 
rijn Franx for useful discussions and for comments on the 
manuscript. 
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APPENDIX A: THE FUNCTION dJx/ds 2 

In order to evaluate c g defined in equation (3.4), we need 
to calculate \dJ\/ds 2 \ at fixed vo and A m . The action J\ 
as a function of the turning points (Vo, Ai, A2) is defined in 



equation (2.11) as a single quadrature. Upon transformation 
to (i/q, A m , e) we have 



A m + e 

j _ _\/2 I dA V(A-A m + e)(A m -A + e) 
A ~ t J (A + a)VA + 7 

x \/ (A-f )^[fo, A m — e, A, A m + e]. 
Since s = e/(A m + a), we have, at fixed A„ 
dJx _ (Xm + af dJ x 



(Al) 



ds 2 



2c 



de 



(A2) 



The integrand in equation (Al) vanishes at the lower and 
upper limits of integration, so we can simply carry out the 
differentiation with respect to e inside the integral. This is 
straightforward upon repeated use of the definition (2.9) of 
divided differences. The result can be written compactly as: 



dJ\ _ (Xm + a) 
ds 2 



A 2 

2 ' dA 



X — vo I f7[fo,Ai,A,A 2 ] 



tta/2 7(A + «)V A + 7 V (A 2 -A)(A-A0 
f (A 2 -A)(A-A 1 )[/[ t / ,A 1 ,A 1 ,A,A 2 ,A 2 ] 'j 

l u\v ,x 1 ,x,x 2 \ y 



This is a function of vq , Ai and A 2 , and hence depends on vo , 
Xm and e, or s. We remark that the expressions for J\ and 
dJx/ds 2 are invariant under the exchange Ai <-> A 2 . This 
means that both these functions are even in s, and hence 
are functions of s 2 . 

We found it convenient to evaluate J\(vo,X m ,s) and 
dJ\/ds 2 (vo, A m , s) by transformation to the integration 
variable w, defined as 



2A — Ai — A 2 A — A^ 



Ai+A 2 



(A4) 



Then dX = s(X m + a)dw, and the integration limits are 
w(Xi) = —1 and w>(A 2 ) = 1. As a result 



J\ = — s 2 (A m + a)\/A-f [dw\/T- 
7r J 



\/l + (l-£o)sw jU[v , Ai, A, A 2 ] 



dJ\ _ (A m + a)\/A^ 
ds 2 ~ 7:^/2 



1 + SM) 
1 



A + 7 



-fo 



dw 



(A5) 



\/l + (l-£o)sw / U[v ,Xi,X, A 2 ] 



1 + 
\2/-| „.. 2 



A + 7 



f s 2 (A m + a) 2 (l-«> 2 Mt/o,Ai,Ai,A,A 2 ,A 2 ] ) 
X l U[v ,X 1 ,X,X 2 ] y 

where we still have to substitute A = A m + sw(X m + a), 
Ai = A m — s(A m + a), and A 2 = A m + s(A m + a). The 
quantity xo is defined in equation (4.19): 



x 



-a — vo 



(A6) 



so that < xo < 1. It is a constant as far as the integration 
over w is concerned. 
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In the thin-orbit limit we have Ai = A m = A2, i.e. 
s = 0, and hence 



Jx = 
dJx 



ds 2 



(A m + a) 



(A m — vo)U[i>o, A m , A m , A m ] 



(A7) 



2(A m + 7 ) 



Substitution in equation (3.4) for c g now immediately gives 
(3.5). 



APPENDIX B: APPROXIMATIONS AT LARGE 
RADII 

The various third-order divided differences U[to, t\ , T2, T3] 
that occur in the fundamental integral equation simplify 
when the potential V becomes Keplerian cx -G¥/A -1 ' 2 
when A — > 00 (ZH, Hunter & de Zeeuw 1992). Here we 
need the case To < — a and —a <C n, T2, T3. Upon substitu- 
tion of the asymptotic behaviour U(X) ~ -GMX 3/2 in the 
definition (2.9) we obtain 



U\(J,V,X\ ,A2] 



GM 



(VAi+ VA 2 ) (VAi+ VA 3 ) (VA 2 + VA 3 ) 
GM 



(Bl) 



t/AiA2(t/AT+a/A 



The function f7*(A, v\ vq, Ai, A2) defined in equation (2.15) 
can therefore be approximated by 



U* 



(GMf 2 (^x- 1 +Vxfi 2 (Vx+VY 2 f 



/2 



(B2) 



4 a/aTMa/aT+a/aT) 7 / 2 

In terms of the variables s and t this can be written as: 
?7* ~ v _„/„ — v — L(s,t), (B3) 



2 9 / 2 



A 9 / 4 



where 



)/2 [a/T+I + a/T^] 1/2 [VT+t + \/T+7\ 



L(s,t) = 2 



1 1/2 



a/1-s 2 [a/T 3 ! + \/T+I 



7/2 



(B4) 



so that i(0, 0) = I and £(s, t) is even in s. 

When —a <C Ai < A2, the integral (2.11) for the action 
J\ is elementary and independent of i/q. It is given by 



, a 1 / 4 \!/ 4 \2 

a/2GM (A2 " Al } 



JX (a/AT+a/AT) 1 / 2 ' 

Transformation to the variables s and t results in 

V /4 [ (l + s)1 / 4 _ (l _ s)1 / 4 ] 2 



Ja ~ a/2GM- 



,1/2 



(B5) 



(B6) 



(1+*) 1 ' 4 [yr+i + yr37 j 

Straightforward differentiation with respect to s now gives 
dJ x \/GM A 1/4 



9s 2 4 (1 + i) 1 /' 

where 

a/2[4 + 2a/1-s 2 



■aw 



(1- s 2 ) 1/4 (a/T^+a/T+^) 



(1-s 2 ) 3 / 4 [• x /r :i s + \/T : n 



1 5 / 2 



(B7) 



(B8) 



so that h is even in s. It is not difficult to show that h(0) = 1 
and h(s) > 1 for < s < 1. 



The normalization function c g (yo, A m ) can now be eval- 
uated by substituting the above approximations in the def- 
inition (3.4). It becomes independent of vq, and can be 
written as 

4 A 5/4 

* „ ^TT„ ,^M - (B9) 



C„a/GM (l+i) 5/4 



phere the constant C g is given by 



C g = ds g srn (s)h(s) 



(BIO) 



Since h(s) > 1, it follows that C g > 1 for normalized </ sm . 

The intrinsic velocity moments are computed by insert- 
ing (3.18) as weight functions in the fundamental equation 
(3.17). With the help of equation (Bl) we can approximate 
the velocities (3.18): 

2 GM \. 

V\ = —J=rL v (s,t), 



V v = —j=rL V v (s,t)(\ - U), 

2 GM 4, 
V<p = —j=rLZ(s, t)U, 

where the (s, ^-dependent part has been separated: 

1j v — „ X 



(Bll) 



1j v ■ 
1j v ■ 



(1 + ^2(71^7+71^) 

1 

(VT+t + VT^)(VTTt + VTT. 

2((1 + t) 2 -s 2 ) 



(B12) 



VT 



sHVT 



' + VT+~s)VTTt' 



2a/T 



(a/i - s + VT+~s~)VTTt' 

so that Lt(0, 0) = L v v (0, 0) = 1 and £*(0, 0) = 0. 

APPENDIX C: APPROXIMATIONS NEAR THE 
FOCAL CORNER 

Near the focal corner in the (vo, A m )-plane, where A m = 
vo = —a, the function U\vq, Ai , A, A2] can be approximated 
by U[—a, — a, — a, — a] = U"'( — a) > 0, so that it can be 
taken out of the integral for Ja. It then follows that 



dJ\ (X m + a)^X m -v r— 

~fi~2 , ^U'»(-a)j(x ,s), (CI) 

a/2(7 — a ) 



phere 



N 1 / dw Vl + {l-x )s 
j(x ,s) =- 



w J VT- 



1 + sw 



(C2) 



The trigonometric substitution w = cos t, followed by use of 
the integral tables of Byrd & Friedman (1971), shows that 

j(x ,s) = 2 [(l_ Xo ) A -(fc)+^n(a 2 , *)] ,(C3) 

wy/l + {l-x )s 1 + s 

where K and n are the complete elliptic integrals of the first 
and third kind, respectively, with arguments given by 



2s 

T+~s' 



2 _ 2(\ — X )s 

1 + (1 — x )s 



(C4) 
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In the thin-orbit limit j(0, 0) = 1, so that expression (CI) 
reduces to (A7) evaluated at A m = vo = —a. Two special 
cases of interest are 



APPENDIX D: THE INTEGRAL (4.22) 

When x = we have xo = 0, and equation (4.22) reduces to 



3(0, s) 

J(M) 



2 



1 



K(k) = 2 F 1 (- 



; 1; s 



(C5) 



Here we have used formulae 8.114 and 9.134.1 of Gradshteyn 
& Ryzhik (1980, hereafter GR), to write j(0, s) explicitly as 
a (hypergeometric) function of s 2 . 

The function c g can now be approximated as 



c g (vo, \ r , 
with 

Jg(x ) = 



\/2(7-«) 



1 



yJU"'{-a) J g( x o)' 



ds 2 g sm (s)j(x ,s). 



(C6) 



(C7) 



1 



The M-integration is a factor, and is easily evaluated as 
7r upon the trigonometric substitution u = cos 2 z. The 
^integral equals 2£(A;)(1 + s)" 1/2 (l - s), where E(k) is 
the complete elliptic integral of the second kind and k 2 = 
2s /{l + s). We express it as a hypergeometric function by 
means of formula GR 3.133.15, and use GR 9.134.1 to write 
it as a function of s 2 . The result is 



p / n x 1 f Ssm(s) 1 j 2 



(D2) 



The case x = 1 can be done in a similar way. Now xo = 1, 
and we have 



This shows that c 9 (i/o,A m ) has radial behaviour near the 
focal corner: its value depends on the direction along which 
the focal corner is approached. The function J g (xo) equals 
1 for all xo when c/ sm = <*>(s 2 ), and it varies slowly with xo 
for sharply peaked </ sm . 

For s max < 1 we can express J g (xo) in terms of the 
moments (3.7) of c/ sm by expanding j(xo, s) in powers of s 2 . 
We give the result for the two special cases x = and x = 1: 



Fg{1) = —Z[y£jt= 7 2 i g , m (s)dt 



The M-integration now gives tt/2, while the t-integral fol- 
lows from GR 3.163.1 after the substitution sx = cos z, and 
results in 7r(l — s 2 ) -1 '' 2 . This leaves 



J «(°) = ^7fL *j < s ^ 

A: = 



■h( 



CO 



(C8) 



2A: \ 
S Jsm). 



J g (0) and J g (l) can be evaluated explicitly for the functions 
c/sm defined in equation (3.6). Use of GR 7.512.11 gives 



J g (0) = 2 J Fi(i,f;2 + g; S Lx) 
J g (X) = 2^1(1,1; 2 + g ;sL) 



(C9) 



These hypergeometric functions equal one for s max = 1, and 
are well-behaved for q > and < s max < 1. For q = we 
can write 



J 9 (0) 
J 9 (l) 



8A/T+1 



37TS 



max 

2 



■[£(*)-(!- w)#(A;)], 



(CIO) 



I + a/T 



where E is the complete elliptic integral of the second kind, 
and k 2 = 2s m ax/(l + Smax). When q = we have 



J 9 (o) 

J 9 (l) 



r(2+g)r(i+ g ) 
r(K?)r(!+?) 
1+1 

l + q' 



These expressions are valid for s m ax 
case J 9 (0) = 16/3tt\/2 and J g (l) = 2. 



(Gil) 



1 as well, in which 



p (r\ - 1 /" ds 2 ffsm(s) 



(D4) 



For Smax < 1 we can express these integrals in terms of 
the moments (s 2k g srn ) of c/ sm defined in equation (3.7) by 
expanding the integrand in powers of s 2 . This gives 



^ 9 (0) 



i 4 ~ r(*+f)r(*+f) 



J 9 (0) tta/2 



E 



A: = 



(S J.m), 



(D5) 



The integrals (D2) and (D4) can be evaluated explicitly for 
the specific choice (3.6) for the function g sm (s). Use of GR 
3.197.3 gives 



W = j}oj 2F ^lh 2 + ^ s ™x 

F g (l) = j-^-y 2^1(1,1; 2 + g; S Lx) 



(D6) 



with J g (0) and J g (l) given in equation (C9). These hy- 
pergeometric functions equal one when s m ax = 0, and are 
well-behaved for q > — 1 and < s m ax < 1. For s m ax = 1 
we obtain 



p rtn i F(2+g)r( g ) 

9() J g (0)r(f+g)r(|+g)' 

F ' il) = m-ir' (g>0) 



(D7) 
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When q = we find 

8[K(k)-(l + s m ^)E(k)] 



tA/l + s max Jg(0) 



— ln(l-Smax), 



(DJ 



with A; 2 = 2s max /(l + s max ). It follows that F g (0) and F g (l) 
diverge logarithmically when s max — ► 1, but their ratio ap- 
proaches 3/2. 
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